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Abstract 



A lattice QCD calculation of the kaon B parameter Bk is carried out with 
the Wilson quark action in the quenched approximation at /? = 6/5^ = 5.9 — 
6.5. The mixing problem of the As = 2 four-quark operators is solved non- 
perturbatively with full use of chiral Ward identities employing four external 
quarks with an equal off-shell momentum in the Landau gauge. This method, 
without invoking any effective theory, enables us to construct the weak four- 
quark operators exhibiting good chiral behavior. Our results for Bk with 
the non-perturbative mixing coefficients show small scaling violation beyond 
the lattice cut-off ~ 2.5GeV. Our estimate concludes 5;^(NDR, 2GeV) = 
0.69(7) at a"^ = 2.7 — 4.3GeV, which agrees with the value obtained with 
the Kogut-Susskind quark action. For comparison we also calculate Bk with 
one-loop perturbative mixing coefficients. While this yields incorrect values 
at finite lattice spacing, a linear extrapolation to the continuum limit as a 
function of a leads to a result consistent with those obtained with the Ward 

identity method. 

11.15.Ha, 12.38.Gc, 13.75.Cs 
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Typeset using REVT^ 



I. INTRODUCTION 



The kaon B parameter defined ratio 

l{K^\s^,^,d\0){0\-s^,^,d\K^) ^ ^ 

is one of tlie fundamental weak matrix elements which have to be determined theoretically 

for deducing the CP violation phase of the Cabibbo-Kobayashi-Maskawa matrix from ex- 
periments. Lattice QCD calculation is expected to evaluate Bk precisely incorporating the 
long-distance effects of QCD. Much efforts have been devoted over the years for this purpose, 
using both the Wilson and the Kogut-Susskind(KS) quark actions. Successful calculations 
of Bk have been achieved so far with the KS quark action, taking advantage of the correct 
chiral behavior of the matrix element ensured by U(l) chiral symmetry [0,0, while studies 
with the Wilson quark action are rather stagnant. There are two purposes for us to try 
to advance the calculations of Bk with the Wilson action. One of them is to verify the 
consistency between the Wilson and the KS results, which would give a full credit to the 
lattice QCD calculation. The other is an application to the heavy-light system for which 
the interpretation of flavor quantum numbers with the KS action is difficult. 

An annoying defect of the Wilson quark action is explicit breaking of chiral symmetry at 
finite lattice spacing. For the calculation of Bk the problem appears as a non-trivial mixing 
of the weak As = 2 four-quark operator of purely left handed chirality with those of mixed 
left-right chirality. Early studies showed that the mixing problem is not adequately treated 
by perturbation theory, leading to an "incorrect answer" for the matrix element 0. Most 
calculations of Bk have then tried to solve the mixing problem non-perturbatively with the 
aid of chiral perturbation theory [^0, and have succeeded in giving reasonable estimates 
for Bk- This method, however, is not promising from a point of view to control systematic 
errors, since it contains large uncertainties from higher order effects of chiral perturbation 
theory which survive even in the continuum limit. 

An essential step toward a precise determination of Bk is to control the operator mixing 
non-perturbatively without resort to any effective theories. The failure of the perturbative 
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approach suggests that higher order corrections in terms of the couphng constant might be 
large in the mixing coefficients. Presence of large corrections in powers of the lattice spacing 
a in the mixing coefficients is also a possibility. In order to deal with this problem, the Rome 
group has proposed the method of non-perturbative renormalization(NPR) 0. Numerical 
results based on this approach show an improvement of the chiral behavior of the As = 2 
operator 0]. 

In this paper we propose an alternative non-perturbative method to solve the operator 
mixing problem which is based on the use of chiral Ward identities [§]. This method fully 
incorporates the chiral properties of the Wilson action explicitly. We also reexamine the 
question if perturbative mixing coefficients lead to erroneous results for in the continuum 
limit. Our simulations have been made within quenched QCD aX {3 = 5.9 — 6.5 keeping the 
physical spatial size approximately constant at 2.4fm. Chief findings of our calculation have 
already been presented in Ref. and we give in this article a detailed description of the 
implementation of our method and the results of our analyses. 

This paper is organized as follows. In Sec. y we describe the formalism of our non- 
perturbative method to determine the mixing coefficients for four-quark operators based on 
the use of chiral Ward identities. The perturbative expressions for the overall renormalization 
factors are also given. In Sec. |TTl| we present our data sets and give a description of the 
calculational procedure for B^- Results for the mixing coefficients are given in Sec. where 
we compare our results with those for the NPR method. The overall renormalization factors 
determined by the NPR method are compared with those obtained by the perturbative 
one in Sec. 0. In Sec. ^ we examine the chiral properties of the four-quark operators 
constructed with the Ward identity method. Final results for Bk are presented in Sec. [VII . 
Through Sees. 0, |V], |VI|, and |V11| we also present results with the perturbative method for 
comparative purposes. Our conclusions are summarized in Sec. |VI11| . 

II. FORMULATION OF THE METHOD 
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A. Determination of the mixing coefficients 



We first derive the generic form of the chiral Ward identities in a standard manner [Q. 
The Wilson quark action is defined by 

+ ^V^(x)(mo + 4)?/'(a;), 



where ip = {u, d, s) represents the up, down and strange quark fields. The conventional 
hopping parameter is given by 

1 



K 



2mo + 8 

Under the flavor SU(3) chiral variation defined by 

- A" 



(3) 



(4) 
(5) 



with A" (a = 1, . . . , 8) the flavor matrices normalized as Tr(A'^A^) = 26"'^, the naive Ward 
identity that follows from the Noether procedure takes the form: 



{V,A^;''^''{x)Oix^,■■■,Xn)) 
2mo{P\x)0{xi, ■ ■ ■ , Xn)) + (X'^(x)C>(xi, ■ ■ ■ , 



(6) 



Xn)) 



where the pseudoscalar density P*^, the extended axial vector current A'f and its divergence 



are defined by 



A" 



AT'\^) = 2 



'^{x)Y^)^^5Uf,{x)^{x + fl)+ ^{x + /i)y7M75f/i(a;)^(a;) 



(7) 
(8) 
(9) 



and the X" term is given by 



A" 

+{x ^ X- fi)] + 8^/'(x)y75^(a;). 
The X" term mixes with and V^^A^ under renormahzation. Thus we write 

X^ix) = - 25moP"(x) - (Z^"*'" - l)V^^^"*'»(x), 

where should satisfy: 

(i) On-shell matrix elements vanish in the continuum limit, 

{a\X''{x)\l3) = 0{a). 

(ii) Off-shell Green functions have only contact terms up to terms of 0{a), 

{X''{x)0{x,, x^)) =Y.^{x- Xi){0','^{x^, + 0{a). 

i 

Defining the renormalized axial vector current by 

and the renormalized quark mass by 

m = mo — Sruo, 

the Ward identity takes the following form: 

{V,Af'''{x)0{x„...,Xr.)) 
= 2m(P«(x)C»(xi, . . . , + {X^{x)0{x,, . . . , Xn)) 

For finite quark masses it is also useful to take a four-dimensional sum over x, which 
2m ^(P«(x)C»(xi, . . . , + Y,{X-{x)0{x^, x„)> 

X X 

+ ^^(TO(a;i,...,x„)) = 0. 



We note that X°(x) in eqs.(0) and (|T^ generates only contact terms up to terms of 0(a). 

Let us consider a set of weak operators in the continuum {Oi} which closes under flavor 
chiral rotations 6°'Oi = ic^^Oj. These operators are given by linear combinations of a set of 
lattice local operators {Oa} as Oi = J2a ZiaOa- We choose the mixing coefficients Zia such 
that the Green functions of {Oi} with quarks in the external states satisfy the chiral Ward 
identity to 0{a). This identity is obtained from eq.(|l7D: 



X k k 

-^T.P^iO)]l^P{pk)S^^P{pl)) + 0(a) = 0, 

I ki^l 



where is the momentum of the external quark. We note that the first term in eq. 
comes from the chiral variation of the Wilson quark action and the third represents the 
chiral rotation of the external fields. Since the identity (|T^) is linear in the renormalization 
constants, the overall renormalization factor cannot be fixed. Furthermore, with quarks in 
the external states, calculations have to be made in some fixed gauge, e.g., the Landau 
gauge. 

The 0(a) term is governed by the typical QCD scale Aqcd at low external quark mo- 
menta, while powers of p^ct become the dominant source of cut-off effects as momenta in- 
crease. To be able to impose the Ward identity to 0(a), we need to restrict the external 
momenta by the condition ^ 1/a for large momenta. On the other hand, no such bounds 
exist for small momenta for the validity of the identity itself as long as Aqcdo ^ 1. 

The parameter = (mg — 5mo)/Z^^* in eq.(|TBp is determined from the PCAC(partial 
conservation of axial vector current) relation obtained from an application of the Ward 
identity for O = P^{y) |0|: 



{V ,Af^\x)P\y)) = 2p^{P\x)P\y)) (19) 

1 



1 ^nh.^ inh^ A 



3 2 



i,{y)) + 0{a). 



For the determination of we employ another Ward identity: 
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where we take O = A^^{y)V^{z) in eq.(|T^ |T^. The local axial vector current and the local 
vector current are defined by 

Al{y) = ^{y)YW^iy)^ (21) 
V:{z)=^{z)^^Mz), (22) 

and Za and Zy are the renormalization factors for v4[^(?/) and Vp{z), respectively. The 
identity ( PD| ) can be regarded as a set of equations for Z'^^ and Za/Zy. Two independent 
equations are obtained from fi = u = 4 and fi = u = i (i = 1, 2, 3). 
The continuum four-quark operator relevant for Bx is given by 

Ovv+AA = {s'j^d){s'j^d) + (s7^75C?)(s7/.75C?) (23) 

where () means color trace, and the parity violating part of the operator which does not 
contribute to is dropped. To fix the mixing coefficients for the lattice four-quark opera- 
tors, we may choose a particular SU(3) fiavor chiral rotation to be applied for Ovv+aa- In 
order to avoid complexities in numerical simulations it is essential to avoid fiavor rotations 
that yield operators which have Penguin contractions and hence mix with lower dimension 
operators. Assuming SU(2) symmetry m„ = we employ the = diag(l,— 1,0) chiral 
rotation, under which Oyv+AA and Ova = (s7^(i) (37^756?) form a minimal closed set of the 
operators: 

6^^dvv+AA = -tOvA. (24) 
d^'OvA = -i\dvv+AA. (25) 

Since Ovv+aa and Ova are dimension six operators with As = 2, we can restrict our- 
selves to dimension six operators for the construction of the lattice operators corresponding 
to them. The set of lattice bare operators with even parity is given by 
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VV = (s7^(i)(s7^(i), 
AA = (s7^75rf)(s7/,75rf), 
SS = {sd){sd), 



and the set with odd parity is 



VA = (s7^(i) (57^756?), 

SP={sd){s-f5d), 

TT = ^{sa^^d){saf,^-f5d), 



(26) 
(27) 
(28) 
(29) 
(30) 

(31) 
(32) 
(33) 



where a^i, = [7^, 7i/]/2. We rearrange these operators into the Fierz eigenbasis, which we find 
convenient when taking fermion contractions for evaluating the Green functions in eq. ([T8|) , 







VV + AA 




(+. 


+), 


(34) 


0,= 




SS + TT + PP 




(+. 


+), 


(35) 


02 = 




SS - \TT + PP 






+), 


(36) 


03 = 


{VV 


-AA) + 2 {SS - 


PP) 




+), 


(37) 


04 = 


{VV 


-AA)-2 {SS - 


PP) 




+), 


(38) 



05= VA (+,+), (39) 
0^ = SP+\Tf (+,-), (40) 
Oj = SP~\Tf (-,-). (41) 

Here the first sign after each equation denotes the Fierz eigenvalue and the second the CPS 
0] eigenvalue. The Fierz eigenbasis we employ is different from that chosen by the Rome 
group 0] based on one-loop perturbation theory. 

The parity odd operators 06,7 are CPS odd while 05 is CPS even, and hence 05 does 
not mix with 06,7 under renormalization, where we assume irtd = iris in the quark action. 
Therefore the mixing structure of these operators is given by 
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Ovv+AA 
Zvv+AA 

Ova 

ZvA 



Ovv+AA = ZqOo + ZiOi + ■ ■ ■ + 2:404, 

Ova = Z5O5, 



(42) 
(43) 



where Zw+aa and Zy^ are overall renormalization factors, and we take zq = 1. 

Let us consider an external state consisting of two s quarks and two d quarks, all having 
an equal momentum p. Under chiral rotation the Ward identity (|18D for such an external 
state takes the following form: 



Fvv+AA = -2p^ZJ*5](P»(x)^6yy+^A(0)s(p)5(p)J(p)J(p)) 



(44) 



-{dvAiOrsip)~sip)dip)dip)) 



-{2Ovv+aa{0)s{p)s{p) 



dip)) 

) + 0{a) = 0, 
FvA = -2p^ZTY.^P\x)dvA{mP)s{p)d{p)d{p)) 



-{-Ovv+AA{0)s{p)s{p)d{p) 



d{p)Y 



(45) 



-{-dvv+AA{0)s{p)s{p)d{p)d{p)) 



-{OvAio)sip):sip) 
-{dvA{ors{p)Hp)d{p) 



~d{p)^ 



d{p) 



dip)) 

75 



) + 0(a) = 0. 



We obtain the amputated Green functions for Fyy+AA and Fva by truncating the external 
quark propagators according to 



VV+AA 



G;\p)G;\p)Fvv+AAGf{p)Gfip), 



TvA = G;\p)G;\p)FyAGj'ip)Gj'ip), 



(46) 
(47) 



where denotes the inverse quark propagator with the flavor q. 

Let Pi {i = 0, ... ,7) he the tree-level Dirac components corresponding to the four-quark 
operators in the Fierz eigenbasis Oi {i = 0, . . . , 7), e.g., 



pafSSX ^ .A ^ (7,75)"^(7.75) 



sx 



(48) 



Since QCD conserves parity one can write 
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= T,P„ (49) 

^VV+AA 

= ToPo + TiPi + ■ ■ ■ + T^P,, (50) 



^VA 



ZvA 

where the c-number coefficients Fq, . . . , Ts are obtained by applying the suitable projection 
operators to Tw+aa/ Zvv+aa and Tva/ Zva, e.g., 



-^PaXS _ 1 

~ 128 



P' 



l^H' + (757.)''"(757.)'1 , (51) 



corresponding to P^^^^ . Expressing Oyv+AAyA in eq.(|T8D in terms of lattice operators, we 
obtain six equations for the five coefficients zi, . . . , ^5: 

Ti = 4 + c\zi + ■■■ + d,Z5 = 0(a), i = 0,...,5. (52) 

This gives an overconstrained set of equations, and we may choose any five equations to 
exactly vanish to solve for zf. the remaining equation should automatically be satisfied to 
0(a). We choose four equations to be those for z = 1, . . . , 4, since Oi, . . . , O4 do not appear 
in the continuum. The choice of the fifth equation, i = or 5, is more arbitrary. We have 
checked that either Fq = or Fs = leads to a consistent result to 0(a) for zi, . . . , Z4 in the 
region pa^l. In the present analysis we choose F5 = 0. 

Let us remark here that the equations obtained in the NPR method corresponds 
to Fj = for i = 1, ... ,4 in which the contributions of the first term due to the quark 
mass contributions and the third term representing chiral rotation of quark fields in the 
Ward identity (p!8[) are dropped. The authors of Refs. P,|7|JTT| argued that the NPR method 
is equivalent to the Ward identities in the limit of large external quark momentum p. The 
reasoning is that the first term is suppressed by one power of p due to the explicit quark mass 
factor, and the third term does not yield chiral-breaking components since the inverse quark 
propagator for large p has the form oc iJ2^j.1^iPfi which anticommutes with 75. Under 
these circumstances the first, third and fourth terms of Fva in eq.(^5D become irrelevant 
for the determination of the mixing coefficients zi, . . . ,Z4. However, the latter point is not 
correct at finite lattice spacing. The inverse quark propagator does not anticommute with 

11 



75 in the large momentum region because the contribution of the Wilson term in the quark 
propagator becomes larger, and hence not negligible, as the momentum increases. Therefore 
the third and fourth terms of Fva in eg. (|45|) yield components having Dirac structures other 
than VV + AA after truncating the external quark propagators. In conclusion the NPR 
method is a part of the Ward identities; at the large external quark momentum p the former 
becomes equivalent to the latter up to 0{pa). In Sec. |V| we show the difference between 
them numerically. 

B. Matching of lattice and continuum operators 

In our earlier report we employed the NPR method of Ref. to evaluate the overall 
renormalization factor Zyy+AA in eq.(^). For the reasons discussed in Sec. |V] we use the 
perturbative estimate in the final analysis presented in this article. 

One-loop perturbative renormalization of the = 2 operator is written in the following 
way ill: 

Ovv+AA = Zvv+aaOo + ^Z* ^-Ci - -O3 - , (53) 

with Oi (i=0,. . . ,4) being the Fierz eigen operators defined in eqs.(|^) — (^). Employing 
the MS subtraction scheme with naive dimensional regularization(NDR) for the continuum 
theory the renormalization factors are given by [0,^ 

4 \n{fia) + Avv+aa) , (54) 
for NDR, (55) 

(56) 

where /i is the renormalization scale. The diagonal part Ayv+AA is affected by the renor- 
malization scheme in the continuum, while the mixing part Z* is independent. We collect 
the value of Aw+aa for the dimensional reduction(DRED) scheme in the Appendix. 

Including the normalization of quark fields v^8!K7y^l72^^^-"378^ HHHHl tadpole- 
improved by the factor uq = {8Kc)~^ leads to 
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Zvv+AA = 1 + 7^ (■ 

Avv+AA = -50.841 
Z* = 9.6431, 



'VV+AA 



(± 



a. 



1 + — (-4 ln(/xa) + Avv+AA + 2n x 5.457) 
47r 



\2K 8K^ 

Here is the critical hopping parameter where the pion mass vanishes. We use 

^ 1 - 5.457a,/4 



(57) 



(58) 



in Ref. for the perturbative estimate of Kc 



With the use of Zyy^j^^ we convert the matrix element on the lattice into that in the 



continuum NDR scheme renormalized at the scale fi = 1/aGeV |ll6|,|5[|: 

5i^(NDR, 1/a) 



Z^.B^,.{K^\Ovv+aa\K') 



WV+AA 
^NDR 



(59) 



where Ovv+aa is defined in eq.(42), and Za is the renormalization factor for the axial vector 



current = 57^756?, which is expressed as pTlJTS 



Za 



-21.061 



for NDR. 



(60) 
(61) 



The value of for the DRED scheme is given in the Appendix. With the tadpole- 
improvement the expression (pO) becomes 



Za 



(. 



1 



a. 



1 + (Aa + TT X 5.457) 
47r 



(62) 



\2K 8K, 

The continuum value at a physical scale fi = 2GeV is obtained via a two-loop renormal 
ization group running from yU = 1/aGeV: 



,(0) 



5^(NDR,/i) 



"ms(1/«). 



2Pl 



(NDR, 1/a), 
(63) 



where /3o,i are the leading and next-to-leading coefficients of the (3 function and 7^'^'^) are 
those of the anomalous dimension for Ow+aa- We take /3o = H, Pi = 102, 7^°^ = 4, and 



(1) 



7 

of Bk. 



-7 appropriate for the zero-flavor case corresponding to our quenched calculation 
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We define another B parameter to investigate the chiral property of the operator 

Ovv+AA- 



5^(NDR, 1/a) = ZVr>r' sZZ:r.Z ' W 



7NDR / 7^0 

"^VV+AA \^ 


Ovv+AA 




7NDR 8 


(0|P|irO)|2 



with Zp the renormalization factor for the pseudoscalar density P = S'-y^d. The continuum 
value of at 2GeV is obtained by running from /i = 1/ aGeV to 2GeV according to the two- 
loop renormalization group. We use 7p^ = —8 and 7p'* = —404/3 0] for the leading and 
next-to-leading coefficients of the anomalous dimension of the pseudoscalar density in the 
zero-flavor case. The one-loop perturbative expression for Zp with the tadpole improvement 
is given by [p!7| , p!5 



Z 



p 







\2K 


skJ I 



ry 

1 + (8 ln(/ia) + Ap + vr X 5.457) 



(65) 



Ap = -30.128 for NDR. (66) 

The value of Ap for the DRED scheme is given in the Appendix. 

The overall renormalization factors Zyy^j^A, Za, and Zp can be alternatively determined 
by the NPR method [^. The NPR method closely follows what is usually done in the 
perturbative renormalization. The vertex corrections are extracted from the amputated 
Green functions for off-shell external quark states with momentum p in the Landau gauge 
according to 

G:'(j9)G:^(j9)(Oyv'+AA(0)5(p)5lp)dlp)J(p))G,^^(p)G,^l(p) = Ao(j9)Po + ■ ■ ■ , (67) 

G:\p){A,iO):s{p)~dip))Gfip) = K,MP,.,. + • • • , (68) 
G:\p){P{mp)d{p))Gf{p) = KiP)P,s + ■■■, (69) 



where P's are the tree- level Dirac components with Pq given in eq. (|48|) and P-^^^y^^-^^ defined 
by 

^"^5 = (7.75)"^ (70) 
= 75"^. (71) 
14 



We should note that the amputated Green functions for the bihnear operators can have 
extra Dirac components besides their tree-level ones, which originate from contribution of 
the higher dimensional operators. The quark wave-function renormalization factor Zg{p) is 
extracted from the quark self energy: 



where the trace is applied for the Dirac and color indices. In terms of the vertex corrections 
and the wave-function renormalization factor one calculates Zyy+AA, Za, and Zp imposing 
the following conditions 

Zvv+aa{p)Z;\p)Ao{p) = 1, (73) 
"1 



ZAip)Z;\p) 



(74) 



Zpip)Z^\p)A,M = ^- (75) 

This renormalization scheme is called the regularization independent (RI) scheme. In this 
scheme the renormalization constants depend on the external state and the gauge. The per- 
turbative values of the renormalization constants Ayy+AA, ^a, and Ap defined in eqs.(^), 



, and (|65|) for the RI scheme are given in the Appendix. 



III. DETAILS OF NUMERICAL SIMULATION 
A. Data sets 

Our calculations are made with the Wilson quark action and the plaquette gauge action 
a.t /3 = 5.9 — 6.5 in quenched QCD. Table | summarizes our run parameters. Gauge config- 
urations are generated with the 5-hit pseudo heat-bath algorithm. At each value of j3 four 
values of the hopping parameter K are adopted such that the physical point for the K me- 
son can be interpolated. The critical hopping parameter is determined by extrapolating 
results for at the four hopping parameters linearly in 1/2K to = 0. We take the 
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down and strange quarks to be degenerate. The value of half the strange quark mass msa/2 
is then estimated from the experimental ratio m^^/mp = 0.648. 

The inverse lattice spacing is determined from the p meson mass rup = 770MeV. The 
physical size of lattice is chosen to be approximately constant at La ^ 2.4fm. To calculate 
the perturbative renormalization factors, we employ the strong coupling constant at the 
scale 1/a in the MS scheme, evaluated by a two-loop renormalization group running starting 
from l/(yf|jg(7r/a) = {TiUp) / gf^^^ + 0.0246 with (Trf/p) the averaged value of the plaquette. 

In order to calculate the mixing coefficients Zi {i = 1, . . . , 5) with the Ward identity 
method and the renormalization factors Zyv+AA, Za, and Zp with the NPR method, the 
latter for purpose of comparison with the perturbative values, we prepare a set of external 
quark momenta p*^*^ = {Px\Py\p'z''iPt'^) — 1? • • • 5 ~ 40). These momenta are chosen 
recursively according to the condition that the [i + l)-th momentum p*^*+^''a is the minimum 
number satisfying 

{p^'^'^af > V ■ (p»a)^ (76) 
^(m)<p(m)<p(m)^ (77) 

for a given value of the increment parameter Sp2 starting with p'^^^a = (0, 0, 0, 2tx /T) where T 
denotes the temporal lattice size. In the case of multiple choices for the i + 1-st momentum 
we take the momentum that has the largest value of pf^^''. The choice of the value of 6p2 
is listed in Table |. We employ the momentum having p^*'^ f» 2GeV among the ^'-^-''s for the 
analysis of B parameters. 

We estimate errors by the single elimination jackknife procedure for all measured quan- 
tities except for the extrapolation to the continuum limit as a function of a. 

B. Calculational procedure 

Our calculations are carried out in three steps. We first calculate m^, mp, pm and 
Z'^^^ using the hadron Green functions. For this purpose quark propagators are solved in 
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the Landau gauge for the point source located at the origin with the periodic boundary 
condition imposed in all four directions. Following eq. (|19|) we can extract the pm parameter 
from the ratio 

^E.^(V4Af'^(x,t)p3(o,0)) 



2E.-(P=^(x,t)P3(0,0)) 



+ {t^T-t + l) 



(78) 



2(0|P3|vr3) 

by fitting a plateau clS cl function of t, where 



< t < T - 1, 



V^AT'\x, t) = AT''\x, t) - AT''\x, t - 1) 



1 cxt , 3 



ext. 3 / 



(79) 



1 r 



A'^^^'^{x, ^) = ^ ""(^5 t)7475f/4(x, t)u{x, t + 1) + u{x, t + l)7475[/4(x, t)u{x, t) (80) 



-{u ^ d)] , 
P^{x,t) = ^ [u{x,t)-f5u{x,t) - {u 



In order to determine we make a zero-momentum projection in y in eq.( 
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where the flavor matrices and A are deflned by 



^0,1,0^ 



0,0,0 
0,0,0 

7ext,3 



^0,0,0^ 



A- 



1,0,0 
0,0,0 



(83) 



At each time slice t we obtain Z'^^'^ and Z^/Zy from the two independent equations cor- 
responding to the choices fi = u = 4 and fi = u = i {i = 1, 2, 3) in eq. ([8^) . In Table ^ we 
summarize the values of m^r, rrip, pm, and Z'^^ for the four hopping parameters at each j3. 

In terms of pm and Z'^^ we determine the mixing coefficients Zi {i = 1, . . . , 5) according 
to the Ward identity (|TB]). The quark Green functions having finite space-time momenta are 
constructed with the point source quark propagators in the Landau gauge. For calculation 
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of the first term in tlie Ward identity ([T8|), we employ the source method [|T9| to insert the 
pseudoscalar density. 

The Bk parameter is extracted from the following ratio of the hadron-three point function 
divided by the two-point functions: 

— > ^Bk{NBR, 1/a) < t < T - 1, 

where operators are defined by 

OKo{x,t) = s{x,t)'y5d{x,t), (85) 
OKo{x,t) = d{x,t)'j^s{x,t), (86) 

4 

dvv+AAix, t) =J2 Zvv+AAZiOi{x, t), (87) 
i=0 

A{x, t) = Zas{x, t)7475C?(^, t), (88) 

with Zyv+AA and Za given in eqs.(|^) and ([621). The contribution of each operator Oj 
= 0, ■ ■ ■ , 4) to -Bi<'(NDR, 1/a) can be measured by the ratio 

j^t /^N ^ Ex,y,z{^Ko{x, T - l)Zvv+AAZiOi{y, t)O^Ko{z, 0)) 

|E^,^(0i?o(f,r-i)i(y,t))E,-.W,'^)(^J^o(i',o)) 

1 (K^\Zvv+AAZiOi\K^) 
L"" ||(0|A|iro)|2 

The sum of R\{t) (i = 0, . . . , 4) is equal to RA(t). The parameter i?^(NDR, 1/a) defined in 
eq.(0) is obtained from the ratio 

Es,yAOKoix, T - l)dvv^AAiy, t)0^o(i', 0)) 

3 E,,yAOK4x, T - i)P(y, t)) Ey'AHy', t)O^Koiz, 0)) 



'-,A^K-{x,T~l)P{y,t))Y.y^^mf.t)0^ 
— > ^E^(NDR,l/a) 0<t<T-l, 

where the renormalized pseudoscalar density is 

P(f, t) = Zps{x, t)-f5d{x, t) (91) 
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with Zp in eq.(|65D. For calculation of the ratios Ra-, and Rp we solve quark propagators 
without gauge fixing employing wall sources placed at the edges of lattice where the Dirichlet 
boundary condition is imposed in the time direction. 

The value of Bk obtained with the Ward identity method depends on the external quark 
momentum p^^^ at which the mixing coefficients are evaluated. To investigate the quark 
mass dependence and a dependence of Bk we employ the averaged value of Bk over the 
five momenta from p^*"^-* to where p^*^ represents the momentum nearest to 2GeV. 

We employ the same procedure for the analysis of Bk- 

IV. RESULTS FOR MIXING COEFFICIENTS 

In Fig. |l| we plot a typical result for the mixing coefficients (i = 1, . . . , 4) as a function 
of the external quark momenta for the case of = 0.15034 at (3 = 6.3. In order to evaluate 
the mixing coefficients we need to choose a specific scale p^*^ that satisfies the condition 
p''*^a ^ 1 to avoid cut-off contaminations. We observe that the mixing coefficients show 
only weak dependence over a wide momentum range 0.02^p^a^^l.O, albeit zi and Z2 have 
large errors in the small momentum region p^a^^O.l. This enables us to evaluate the mixing 
coefficients with small uncertainties from the choice of the momentum p^*\ We adopt the 
value p^*^ ~ 2GeV, which we find to always fall within the range of a plateau for our runs 
at /3 = 5.9 -6.5. 

Let us compare the mixing coefficients obtained by the Ward identity (WI) method with 
those by the NPR. Since the NPR method does not employ the full Ward identity of eq.(0), 
it is important to investigate differences in the mixing coefficients between the NPR and 
the Ward identity methods. In Fig. ^ we present the result for the mixing coefficients Zi 
obtained with the NPR method. The NPR result shows a strong scale dependence in the 
region p^a^^O.3, which contrasts to the Ward identity result in Fig. |1|. We suspect that this 
behavior of the NPR result originates from physical non-perturbative contributions, which 
survive even in the continuum limit (see also Sec. |V|). Although we observe a similar scale 
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dependence for the two results beyond the scale p^a? ~ 0.3, the numerical values for Z2 and 
^3 show clear deviations beyond the error bars between the WI result and that from the 
NPR for momenta as large as p^a? ~ 2. This is contrary to the expectation that the NPR 
would be equivalent to the Ward identities in the limit of large external quark momenta 

mm- 

The difference between the NPR and the WI methods comes from the first and third 



terms in eg . ([T8|) . To investigate the contribution of each term to the mixing coefficients we 
reevaluate the mixing coefficients using the Ward identities without the first term or the third 
term. The former result is plotted in Fig. |^(a) and the latter one in Fig. ^(b). Comparison 
between Fig. |^(a) and Fig. |l] demonstrates that the contribution of the first term to the 
mixing coefficients is remarkable in the lower momentum region p^a^^O.3. Above this scale, 
the first term seems to play a minor role on the determination of the mixing coefficients. 
On the other hand, comparing Fig. ^(b) with Fig. |l| the essential contribution of the third 
term to the mixing coefficients is observed over a wide range of external momentum, even 
up to p^a? ~ 2. In the Ward identity method we can neglect neither the first term nor the 
third one. 

Figure ^ shows the quark mass dependence of the mixing coefficients zi {i = 1, ... ,4) 
evaluated at the scale p^*'^ (filled symbols) for the case of /? = 6.3. For comparison we also 



plot the perturbative (PT) estimate for Zi (open symbols), which are given in eq . (|53D , i.e., 

(92) 
(93) 
(94) 
(95) 

where ajj^{l/a) is used for the strong coupling constant. We observe little quark mass 
dependence for the mixing coefficients. 

In Fig. ^ we present the a dependence of the mixing coefficients Zi {i = 1, . . . ,4) evaluated 
at the scale p'^*^ employing the heaviest quark mass at each p. We observe that the a 
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dependence of the mixing coefficients determined by the Ward identities is steeper compared 
to that of the PT estimates to one-loop order. The magnitude of each mixing coefficient for 
the WI method varies nearly in proportion to a, which reduces by 50% between nipa ^ 0.4 
and irtpa ~ 0.2. A possible source of this a dependence of the mixing coefficients is the 0{a) 
term in eg. (P^ : contributions of the 0(a) term are absorbed in the mixing coefficients to 
satisfy the continuum Ward identities at finite lattice spacing. 

Comparing the mixing coefficients for the Ward identity method and those of pertur- 
bation theory in Figs. ^ and ^, we note that a large value of Z2 determined by the Ward 
identities sharply contrasts with the one- loop perturbative result Z2 = 0. The magnitude 
of this discrepancy appears larger than that possibly explained by two-loop contributions; 
squaring a typical magnitude of one-loop terms in Fig. ^only yields a value of order ~ 0.001. 
Discrepancies also exist for the other coefficients, albeit less conspicuous in that the PT re- 
sults agree with those of WI in sign and rough orders of magnitude. In particular the 
magnitude of is larger than that for for all values of /3, which is contrary to the 
perturbative result. 

For our study of the B parameter the mixing coefficient z^ for the parity-odd operator 
Ova is not directly relevant. However, it is instructive to examine the scale dependence of z^, 
because it would take the value 2:5 = 1 in the absence of cut-off dependent chiral symmetry 
breaking effects. In Fig. ^ we plot a typical result for z^. We find a scale dependence 
stronger than those oi Zi {i = 1, . . . , 4) for parity-even operators toward large momenta; the 
value of Zr, significantly deviates from unity as the momentum increases, which measures the 
magnitude of cut-off effects. The quark mass dependence of z^ evaluated at p^*^ is shown in 
Fig. 1^. The value of z^ slightly increases as the quark mass decreases. We do not consider 
the strong scale dependence of z^ to be particularly alarming since z^ evaluated at a fixed 
physical scale p^*^ approaches unity toward the continuum limit as shown in Fig. |^. 
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V. RESULTS FOR OVERALL RENORMALIZATION FACTORS WITH THE 

NPR METHOD 

The NPR method is a possible way to estimate the overall renormalization factors 
Zyy+AAi Z a, and Zp. In Fig. |^ we plot Zyv+AAi and Zp in the RI scheme as a 
function of p^a? for the case oi K = 0.15034 at /3 = 6.3. For comparison we also draw the 
tadpole- improved one-loop perturbative estimates (solid lines) in the RI scheme. The NPR 
result for Zyv+AA in Fig. ^(a) shows an agreement with the perturbative estimate in the re- 
gion 0.1^j9^a^^0.5. The dotted curve in Fig. ^(a) represents the Oq contribution to Zyv+AA) 
which is obtained by neglecting contributions of the mixed operators Oi (i = 1, . . . , 4). We 
observe that the contributions of the mixed operators, leading of which is the two-loop 
radiative corrections, are quite small. 

Figure |](b) shows that Za has a strong p^a? dependence below p^a^ ~ 0.3. This behavior 
is contrary to the expectation that Za should be independent of p^a? as the axial vector 
current has no anomalous dimension. The unexpected behavior may be ascribed to non- 
perturbative contaminations due to the pion pole which could give an important contribution 
at low external quark momentum [0. In Fig. ^(c) we observe a large deviation between the 
NPR result for Zp and that from perturbation theory below p^a^ ~ 1. We suspect that 
this discrepancy is also due to the non-perturbative effects from the pion pole. It should be 
remarked that the contamination due to the pion pole is not a lattice artifact, and hence 
survives even after taking the continuum limit. 

We note that our NPR results for Za and Zp are consistent with those of recent studies 
|p0|^2|. In particular evidence for the existence of pion pole contribution in Zp has been 
reported in Ref. for the Kogut-Susskind action and in Ref. I^T] for the Wilson case. 

For the calculation of B parameters the ratios Zyv+aa/ Z\ and Zyy+^^/Zp are more 
relevant. We show the scale dependence of Zyv+aa/ Z\ in Fig. |10|(a) and that of Zyv+aa/ Z"^ 
in Fig. 0(b). Solid lines are tadpole-improved one- loop results in the RI scheme. We observe 
that the NPR result for Zw+aa/ Z'a has a p^a? dependence opposite to that expected from 
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the perturbative estimate. For the ratio Zvy+A^Z-Z^p the NPR result diverges toward lower 
momentum. 

Application of the NPR method requires the existence of a region Aqcd <^ p ^ l/a where 
we can keep under control both non-perturbative contaminations and cut-off effects, which 
is called a "window" in Ref. We find for Zyy^AA-, Za-, Zp, and their ratios that the lower 
bound of the window depends strongly on each operator; we observe p^a^ ~ 0.1 for Zyv+AA 
and p^a? ~ 0.3 for Z^, while it is difficult to find the lower bound for Zp, Zw+aaI 
and Zvv+aa/ Zl,. It is not clear to what extent non-perturbative contaminations can be 
separated out quantitatively. For these reasons we employ the perturbative estimates for 
Zvv+AA, Za, and Zp, rather than those of the NPR method, to obtain the B parameter in 
our final analysis presented in this article. Numerical values of the parameter are little 
affected by this change since the the ratio Zyy+aa/ Z\ has a similar value at p*^*^ among the 
two methods. 

For completeness let us examine the quark mass dependence and the a dependence of 
Zyyj^AAi Za and Zp taking their results at p*^*^ (vertical lines in Fig. ^. In Fig. [TT| we plot the 
NPR results for Zyv+AA^ Za, and Zp together with tadpole- improved perturbative values 
as a function of = (1/ K — 1/ for the case of /3 = 6.3. While we observe little quark 
mass dependence for Zyv+AA and Za, Zp clearly decreases as the quark mass decreases. 
Figure |12| shows the a dependence of Zyyj^AA-, Za-, and Zp evaluated at p*^*) employing the 
heaviest quark mass at each /3. The NPR results for Zyyj^AA and Za are consistent with 
the perturbative ones at four (3 values, while for Zp we observe a large deviation at each (3. 



VI. CHIRAL BEHAVIOR 

Let us examine the chiral property of the operator Oyy+AA using the matrix element 
B^ defined in eq. (|6^, which vanishes in the chiral limit in the presence of chiral symmetry. 
In Fig. |1^ we show the chiral behavior of i?^(NDR, 1/a) for the case oi (3 = 6.3. Numerical 
values of i?^(NDR, 1/a) for the four hopping parameters at each (3 are summarized in 
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Table |T|. The solid lines represent quadratic extrapolations of the WI and PT results in the 
bare quark mass m^a = {1/ K ~ \jK^j2. The extrapolated value at = is consistent 
with zero, demonstrating a significant improvement of the chiral behavior compared to the 
perturbative result plotted with triangles. 

Figure |T1| shows the jp'd} dependence of -B^(NDR, 1/a) extrapolated to vtiq = 0, which 
are evaluated at each scale p^^\ Two solid lines represent the upper and lower bounds of one 
standard deviation error for the PT result. The values of -B^(NDR, 1/a) for the WI method 
are consistent with zero within error bars in the momentum range p^a^^l, albeit the errors 
become larger toward the small momentum region. 

We plot in Fig. |15| the values of -B^(NDR, 2GeV) in the chiral limit as a function of 
lattice spacing. The numerical values are given in Table The result for the WI method 
becomes consistent with zero at the lattice spacing mpa^0.3(a^0.08fm). In the perturbative 
approach with the one-loop mixing coefficients, chiral breaking effects are expected to appear 
primarily as terms of 0{a) and secondly as terms of O (a ■ (7^ (1/a)) and (^((^^(l/a)) for the 
Wilson quark action. A roughly linear behavior of our results for the perturbative method 
indicate the leading contribution from an 0(a) term. Making a linear extrapolation to the 
continuum limit a ^ 0, we observe that the chiral behavior is recovered. This may suggest 
that the 0{a-g'^{l/a)) and 0(5''^ (1/a)) terms in the mixing coefficients left out in the one- loop 
treatment are small or accidentally canceled. 

VII. RESULTS FOR Bk 



We now turn to the calculation of -B/^(NDR, 2GeV). In Fig. |T6| we present the ratio -Ra(^) 
defined in eq.(|^) using the mixing coefficients determined from the Ward identities with the 
external quark momentum p^*^ for the heaviest quark mass(i^' = 0.15034) and the lightest 
one(i^' = 0.15131) at (3 = 6.3. A good plateau is observed in the range 20^t^75. We make 
a global fit of the ratio -Ryi(t) to a constant over 32 < t < 63 for this data set. The three 
horizontal lines denote the central value of i?A'(NDR, 1/a) and a one standard deviation 
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error band. We note that the error of the fitted result is roughly equal in magnitude to 
those of the ratio over the fitted range, while we would usually expect a smaller error for 
the fitted result. This is because the error of the ratio RA{t) is governed by those of the 
mixing coefficients Zi {i = 1, . . . ,4). Numerical values of i?x(NDR, 1/a) for the four hopping 
parameters at each f] are listed in Table [m| . 

For comparison we also show the results for RA{t) obtained with the perturbative mixing 



coefficients in Fig. |17[ We observe a plateau in the range 30^t^65, which is slightly narrower 



compared to the WI results. A global fit of RA{t) to a constant choosing the same fitting 



range as for the WI case yields the value of (NDR, 1/a) given in Table |lTT[ Let us note 
that the PT results have quite small errors compared to those of the WI method. This is 
because definite values are taken for the mixing coefficients in the PT method. 

In Fig. |1^ a representative result for the contribution of each operator Oi (i = 0, ■ ■ ■ , 4) 
to -Bk(NDR, 1/a) is shown as a function of the external quark momentum for the case of 
K = 0.15034 at /5 = 6.3, which is obtained by fitting the ratio -R^(t) {i = 0,---,4) of 
eq.(|5PD with a constant over the same fitting range as for RA{t)- The contributions of mixed 



operators are nearly independent of the external quark momentum in the range p^a^^l.O. 
An important observation is that the value of -B/^(NDR, 1/a) results from large cancellations 
between the amplitudes of the mixing operators ZiOi (i = 1, . . . , 4), each having a magnitude 
comparable to or larger than that of Oq. This is the essential reason why calculations of 
with the Wilson quark action is difficult; the mixing coefficients have to be known accurately 
including higher order effects both in the coupling constant and the lattice spacing. 



We show the quark mass dependence of i?7^(NDR, 1/a) for (3 = 6.3 in Fig. |1^. We observe 
that the results for the PT method seem to diverge toward the chiral limit, while those for 
the WI method stay finite. We expect a different quark mass dependence for the WI and 
PT results: 

5x(NDR, 1/a) = + B^^mg + C^^m^ for WI, (96) 

Sx(NDR, 1/a) = + B^^ + C^'^rrig for PT, (97) 

nig 
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where A, B, and C are unknown constants. These functional forms are based on the following 
assumption for the chiral behavior of the matrix elements near the chiral limit: 

{K''\6vv+aa\K'') oc m, for WI, (98) 
{K°\dvv+AA\K^) oc const for PT, (99) 
(0|4|ir°) oc (100) 

For the WI and PT methods we interpolate the data at the four hopping parameters with 
the forms (|96|) and (|97D , respectively, to ms/2 and obtain the value of -Bx(NDR, 1/a) at the 
physical point. 



We summarize our final results for i?x(NDR, 2GeV) in Table whose a dependence is 



illustrated in Fig. |2y. The method based on the Ward identity gives a value well convergent 
from a lattice spacing of m^a ~ 0.3. Unfortunately the large errors do not allow us to take a 
linear extrapolation to the continuum limit. We may instead take a constant fit of the three 
results at smaller lattice spacings(a"^ = 2.7 — 4.3GeV) and find i?x(NDR, 2GeV) = 0.68(7), 
which is our best estimate for the WI method. 

Since the origin of the large error is traced to that of the mixing coefficients, we attempt 
to develop an alternative method, in which the denominator of (BUI) is estimated with the 



vacuum saturation of the operator Ovv+aa constructed by the WI method: 



where in terms of eq. (|4^) the vacuum saturation of Ovv+aa is rewritten as 

{K'>\Ovv+AA\K')ys = E^«(^V«l^°)vs, (102) 

i=0 

with 

{K'\Oo\K')ys = 1{K'\A\0){0\A\K'), (103) 
{K'\0,\K')ys = 1{K'\P\0){0\P\K'), (104) 
(i?°|02|i^°)vs = Uk'\P\0){0\P\K'>), (105) 



3 
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{K'\Os\K')ys = -^{K'\A\0){0\A\K') - ^{K'\P\0){0\P\K'), (106) 
{K'\0,\K')ys = -1{K'\A\0){0\A\K') + ^{K'\P\0){0\P\K'). (107) 

We refer to this as the Wlys method, with which the fluctuations in the numerator are 
expected to largely cancel against those in the denominator. In fact, errors are substantially 



reduced with the Wlys method as apparent in Fig. The cost is that the correct chiral be- 
havior of the denominator is not respected at a finite lattice spacing due to the contributions 
of the pseudoscalar matrix element. This contribution brings the Wlys result to disagree 
with WI at a finite lattice spacing, but the discrepancy should vanish in the continuum limit. 
A linear extrapolation of the Wlys results in a yields 5;^(NDR, 2GeV) = 0.562(66). 

This linear extrapolation, however, involves a systematic uncertainty arising from the 
chiral symmetry breaking term cp| (0|P|i^'°) p in the denominator, where Cp = 8/3zi+A/3z2— 
8/3^3 + 16/3^4 from eqs.( |104| ) — ( |107| ). Perturbative contributions to cp starts at two- loop 
order of 0{g^{l/a)) as can be checked from the one-loop expressions for Zi (i = 1, . . . , 4) in 
eg. (|5BD . Since the matrix element {0\P\K'^) diverges in proportion to [(7^(l/a)]~^/^^ due to the 
anomalous dimension of the pseudoscalar operator P, cp\{0\P\K^)\'^ receives contributions 
of form [5f^(l/a)]^^/^^ which diminishes only as a fractional power of 1/loga. To assess the 
systematic error associated with this effect, we estimate the two-loop contribution to cp by 
squaring the typical magnitude of the one-loop terms in zf. e.g., \z°"''^~''°°^{ajY^{l/ a))\^0.08 



at P = 5.9 from Fig. ^ We also estimate 

{0\P\K^) rriK 



{0\A\K°) md + nis 

from the PCAC relation 0, which yields cp|(0|P|ir°)|V(8/3)/|(0|A|irO)|2^0.4. Since 



a-jjg{l/ aY^/^^ decreases by 30% between /? = 5.9 — 6.5, over which a decreases by a factor 
2, this fraction should reduce to ~ 0.16 after taking the continuum limit. Taking ac- 
count of uncertainties in the choice of coupling constant and the mixing coefficients at the 
two-loop level, we estimate the chiral symmetry breaking contribution of the pseudoscalar 
density that survives after a continuum extrapolation linear in a to be ^20%. We conclude 
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5A'(NDR,2GeV) = 0.56(7)(11) for the Wlys method. 

An interesting point in Fig. ^ is that the perturbative calculation (PT), which gives 
a completely "wrong value" at a 7^ 0, yields the correct result for Bk, when extrapolated 
to the continuum limit a = 0. This is a long extrapolation from negative to positive, but 
the linearly extrapolated value -Bx(NDR, 2GeV)=0. 622(69) is reasonable compared with 
those obtained with the WI or Wlys method. This linear extrapolation is a credible choice 
because the chiral behavior of the matrix element {K^\Ovv+aa\K^) is linearly recovered as 



we saw in Fig. We have to make a reservation, however, that this long extrapolation 
may bring an error larger than quoted in the extrapolated value due to systematic effects of 
0{a- g'^{l/a)) and 0{g'^{l/a)). The estimation of these systematic errors is too complicated 
because the matrix elements of the mixing operators have quite different absolute values. 

Each of results from the above three methods suffers from statistical and systematic 
errors of 10 — 20% which are comparable in magnitude. Although the Wlys and the PT 
methods have the advantage of small statistical errors, we recognize that this is offset by 
the difficulty to control large systematic errors when attempting a continuum extrapolation. 
We thus conservatively take the resuh of the WI method 5;^(NDR, 2GeV) = 0.68(7) at 
= 2.7 — 4.3GeV as our final estimate of the present work. 



VIII. CONCLUSIONS 

In this paper we have presented a full account of our method based on chiral Ward 
identities to non-perturbatively determine the mixing coefficients of the As = 2 operator for 
the Wilson quark action in lattice QCD. Implementing the method in a quenched calculation 
carried out at four values of lattice spacing, we have demonstrated the effectiveness of the 
method for constructing the As = 2 operator with the correct chiral property. Our final 
result for i?;^(NDR, 2GeV) = 0.68(7) at a^^ = 2.7 — 4.3GeV shows a reasonable consistency 
with i?ii:(NDR, 2GeV) = 0.628(42) in the continuum limit recently obtained with the Kogut- 
Susskind quark action by us 0. 
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The error of our Wilson result for Bk, however, is still too large to convincingly demon- 
strate that the Wilson and Kogut-Susskind quark actions yields the same value in the con- 
tinuum limit. We emphasize that this large error is not due to an intrinsic defect of the 
Ward identity method. It stems from that of large statistical errors of the mixing coeffi- 
cients, which in turn originates from our use of point source in evaluating relevant quark 
Green functions. Recent work shows that a variant wall source method with the momentum 



source for the off-shell quark propagator P2| , |24| would be effective to diminish the errors of 
the mixing coefficients. 

Another technical point concerns the issue of Gribov copies in the Landau gauge. While 
an earlier study p5[ suggests that ambiguities in the choice of the Gribov copies induce only 



small uncertainties comparable to typical statistical errors in current numerical simulations, 
exploring gauge invariant implementation of the Ward identity methods, either employing 
the external hadron states or the Schrodinger functional, which is free from this problem, 
would be worthwhile. 

In recent calculations of using the 0(a)-improved quark action the chiral property of 
Ovv+AA constructed with one-loop mixing coefficients shows much improvement compared 
to the Wilson quark case [^. This observation can be expected on the ground of our 



perturbative results in the Wilson quark action which suggest in Sees. ^ and [VI I| that the 
leading contribution to chiral breaking effects is 0{a). Toward a precise determination of 
Bx the improvement of the quark action is an essential ingredient. 

A very important physics issue is the effect of quenching. With the KS quark action it 
has been observed that the error due to quenched approximation is small Whether 
this is supported by calculations with Wilson action we must defer to future studies. It 
is straightforward to apply our method once configurations are generated with dynamical 
quarks. 

Finally the application of our method for calculations of Bb would be a worthwhile 
attempt since previous calculations of Bb have relied on the mixing coefficients which were 
calculated perturbatively in the massless limit with tadpole improvement. 

29 



ACKNOWLEDGMENTS 



This work is supported by the Supercomputer Project No.32(FY1998) of High Energy 
Accelerator Research Organization(KEK), and also in part by the Grants-in-Aid of the 
Ministry of Education (Nos. 08640404, 09304029, 10640246, 10640248, 10740107). Y.K. is 
supported by Japan Society for Promotion of Science. 

APPENDIX 

The continuum renormalization scheme dependence of the renormalization constants 
A vv+AA, A^, and Ap defined in eqs.(p4[), ( pOD and (p5| ) have been computed for the NDR, 
DRED, and RI schemes by a variety of authors [p|,p8| at the one-loop level. We consider 
it useful to reproduce them in this Appendix using off-shell external quark states in the 
general covariant gauge for quark self-energy and vertex functions. 

We consider the following operators: 

O,.^,. = \ ^''^'^ (^i7.(l - 75)^^) - iM) , (109) 

O,, = 5,, {^{^,4) , (111) 

where 7^ 7^ represents 7^(1 — 75) ® 7/^(1 — 75) and i, j, k, and / label the color indices. 
We note that the operator Ovv+aa defined in eq. (P^D is parity conserving part of O^l^^l 
with ipi = = s and '?/'2 = = d. 

We draw the relevant one-loop diagrams in Fig. (a) the quark self energy, (b) the 
one- loop vertex correction for the quark bilinear operators, and (c) — (h) the six types of 
the one-loop vertex corrections for the four-quark operator. The off-shell momentum for the 
external quark state is denoted as p. The gauge dependence is parameterized by A expressing 
the gluon propagator as b^yjy? — (1 — X)k^kv/k'^. 

Up to the one-loop level the inverse quark propagator and the vertex functions for P = 
7/^' ® 7/^? 7m75' 75 written in the following form: 
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G-^b)=¥-^S«(p), (112) 



Ar(p) = r + ^A«(p), (113) 

where the superscript (i) refers to the z-th loop level. In Table we compile the results 
for S*-^-* and Ap'* obtained by employing the NDR scheme and the DRED one ||2^. The 



reduced space-time dimension D is parameterized by e as D = 4 — 2e, e > 0. We should 
note that the one-loop vertex corrections yield the extra evanescent operators which vanish 
in D = 4 both for the NDR and the DRED schemes. It is meaningless to give results 
without mentioning the definition of evanescent operators, because the constant terms at 
the one-loop level depend on the definition of the evanescent operators. Our choice is as 
follows: 

1 [2- DY 

^ifmj, = 47p757a.(1 - 75) ® 7a.(1 - 75)7np - ^ 7^.(1 - 75) ® 7m(1 - 75), (114) 

<T = ^".^7.75 - ^7^75, (115) 

D 

= V7^ri - 75j ® 7/.U - 75j - 



E^BR = 5^,7.(1 - 75) ® 7m(1 - 75) - -rl,{l - 75) ® - 75), (116) 



where 6^j_i, is the D-dimensional metric tensor. 

From the results for Tj'^'^\p) and J^y\p) (r = 7^^ ® 7^, 7;i75, 75) we can extract the scheme 
dependence of the renormalization constants Ayy+AAi and Ap, which are summarized 
in Table [V^. For the RI scheme, for which the renormalization constants depend on the 
external states and the gauge, we employ the off-shell external quark state with momentum 
= jj? and the Landau gauge fixing. 
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FIGURES 

FIG. 1. Mixing coefficients zi,...,Z4 obtained with the Ward identity method plotted as a 
function of external momentum squared p'^a^ for K = 0.15034 at (3 = 6.3. Vertical line corresponds 
to « 2 GeV. 

FIG. 2. Mixing coefficients zi, . . . ,Z4 obtained with the NPR method. Parameters are the same 
as in Fig. |l|. 

FIG. 3. Mixing coefficients zi, . . . ,Z4 obtained with the Ward identity (|l8|) neglecting (a) the 
first term or (b) the third one. Parameters are the same as in Fig. ||. 

FIG. 4. Quark mass dependence of mixing coefficients zi,...,Z4 evaluated at p*-*-* ~ 2 GeV 
using the Ward identity (WI; filled symbols) method at /? = 6.3. Perturbative (PT; open symbols) 
results are also plotted for comparison. 

FIG. 5. Comparison of mixing coefficients zi, . . . ,Z4 evaluated at p^*^ ~ 2 GeV using the Ward 
identity (WI; filled symbols) method and the perturbative (PT; open symbols) one as a function 
of rupa. 

FIG. 6. Same as Fig. |^ for mixing coefficient z^. 

FIG. 7. Quark mass dependence of mixing coefficient Z5 evaluated at p^*^ ~ 2 GeV using the 
Ward identity method at /? = 6.3. 

FIG. 8. Lattice spacing dependence of mixing coefficient z^ evaluated at p^*^ ~ 2 GeV using 
the Ward identity method. 

FIG. 9. Renormalization factors {a,)Zvv+AA, (k>)ZA, and {c)Zp in the RI scheme obtained with 
the NPR method as a function of external momentum squared p'^a^ for K = 0.15034 at /? = 6.3. 
Transverse lines denote tadpole-improved perturbative estimates. Dotted curve in (a) represents 
the Oo contribution to Zyv+AA- Vertical lines correspond to p^*^ ^ 2 GeV. 
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FIG. 10. Ratios of renormalization factors (a)Zyy_(_yi^/Z^ and (b)Zyy+^^/Zp obtained with 
the NPR method as a function of external momentum squared p^a^ for K = 0.15034 at /3 = 6.3. 
Transverse lines denote tadpole-improved perturbative estimates. Vertical lines correspond to 

PS 2 GeV. 

FIG. 11. Quark mass dependence of renormalization factors at p*-*-* ~ 2 GeV using the NPR 
(solid symbols) method at /3 = 6.3. Open symbols represent tadpole-improved perturbative esti- 
mates. 

FIG. 12. Comparison of renormalization factors evaluated at p^*^ ~ 2 GeV using the NPR 
(NPR; solid symbols) method and the tadpole-improved perturbation theory (PT; open symbols) 
as a function of mpO. 

FIG. 13. Test of the chiral behavior of 5^(NDR, 1/a) for the Ward identity (WI) and pertur- 
bative (PT) methods at /? = 6.3. Solid curves are quadratic extrapolations to the chiral limit. 

FIG. 14. Dependence of -B^(NDR, 1/a) in the chiral limit on the external momentum squared 
p^a^ at P = 6.3. Horizontal solid lines represent the upper and lower bound of one standard devi- 
ation error for tadpole-improved perturbative result in the chiral limit. Vertical line corresponds 
to w 2 GeV. 

FIG. 15. S^(NDR,2Gev) at = for the Ward identity (WI) and perturbative (PT) meth- 
ods as a function of a. Solid line is a linear extrapolation of the perturbative results to the 
continuum limit. 

FIG. 16. Ratio i?yi(t) using mixing coefficients determined by the Ward identity method for 
{a)K = 0.15034 and {h)K = 0.15131 at /3 = 6.3. Solid lines denote the fitted result and a one 
standard deviation error band. 



FIG. 17. Same as Fig. If: for the perturbative method. 
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FIG. 18. Contributions of the operators Oi {i = 0, ... ,4) to Sx(NDR, 1/a) with Zi determined 
by the Ward identity method for K = 0.15034 at /3 = 6.3. Vertical line corresponds to p^*^ ^ 2 
GeV. 

FIG. 19. Quark mass dependence of ^^'(NDR, 1/a) for the Ward identity (WI) method and 
the tadpole-improved perturbative (PT) one at /3 = 6.3. Open symbols are interpolations of data 
to ms/2. 

FIG. 20. S;<r(NDR, 2GeV) plotted as a function of rupa for the WI, Wlys and PT methods. 
Solid lines show linear extrapolations to the continuum limit. 

FIG. 21. One-loop diagrams for (a) quark self energy, (b) vertex correction for the quark 
bilinear operator, and (c)— (h) vertex corrections for the four-quark operator, p denotes a off-shell 
momentum for the external quark state, and i, j, k, and I label color indices. 
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TABLES 



TABLE I. Parameters of our simulations. See text for details. 



(3 


5.9 


6.1 


6.3 


6.5 


xT 


X 64 


32^ X 64 


40^ X 96 


48^ X 96 


No. conf. 


300 


100 


50 


24 


thermalization 


22000 


32000 


45000 


72000 


interval 


2000 


2000 


5000 


8000 


K 


0.15862 


0.15428 


0.15131 


0.14925 




0.15785 


0.15381 


0.15098 


0.14901 




0.15708 


0.15333 


0.15066 


0.14877 




0.15632 


0.15287 


0.15034 


0.14853 




0.15986(3) 


0.15502(2) 


0.15182(2) 


0.14946(3) 


msa/2 = mda/2 


0.0294(14) 


0.0198(16) 


0.0144(17) 


0.0107(16) 


a-i[GeV] 


1.95(5) 


2.65(11) 


3.41(20) 


4.30(29) 


La [fm] 


2.4 


2.4 


2.3 


2.2 


(Tr[/p) 


0.582 


0.604 


0.622 


0.638 


0'ms(1/") 


0.1922 


0.1739 


0.1596 


0.1480 




1.11 


1.11 


1.15 


1.12 




0.9595 


0.5012 


0.2988 


0.2056 


fitting range for m7r, nip, pm, 


12-20 


14-24 


17-27 


20-30 


fitting range for Bk, 


18-45 


24-39 


32-63 


35-60 
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TABLE II. Meson masses, parameter, and renormalization factor for the extended axial 

vector cuiTciil Ml 3 = 5.9 — G.5 in quenched QCD. 



R 

H 


K 




TYl 
llbp 


Tn I Tti 
llh-j^ 1 llhp 






O.iJ 






n AA'^(^ (\\ 


U.OOUl t 


u.uoou i 1 oo 1 






u. ±o ( oo 






n fi497('fi7'l 


U.UO^OOl O.^ ) 


n Q41 ^^ q^i 




n 1 117(18 
U. ±tj / uo 


n '\( 1 9^1 


u.^oy 1 I oo I 


n 71 79('4'i'l 

U. / J- / ^O J 


n n7fii 7i'\^ 

U.U 1 UJ- i I OJ- I 


U.c/^Ol ±0 J 




u. ±ouo^ 




n ^181 ('9^'! 


n 7fi87('^^'l 
u. 1 uo 1 1 oo J 


U.Uc/OOUl ^i? 1 


U.C/±i7l J-U J 


fi 1 

u. ± 






n ^9^M 9"! 




fl n99^Q('9Q'\ 


U.J7 1 Ul OU J 








U.O^U 1 I UU I 


n fii fiM 1 


n n'^7'^9('9Q'l 

U.UO / O^l ^c/ I 


U.c/OOl ^O i 




U. lOOOO 




U.OUOOl 1 


U.UOOOl Ol ) 


n nii97fi('9Q'\ 


u.yooi J- i j 




1 "^9^7 

U. ±0Z(0 1 




U.OOc/^l OJ- ( 


u. 1 oool uo J 


n nfi778('9Q'l 


n c^'^'^(^A\ 

U.iyOOl J-^ J 


LI. O 


n 1 ^1 SI 


n 1 989('93'l 


n '?'^Al^A^ 


^n4('97'l 


n 01 79^("^4'l 


Q81 




U. ±OUc/0 




u.^u^oi uo J 


R91 M ^ 








U. lOUUU 






RQl M 1 ^ 

U.uy 1 J- -Ly 








0.15034 


0.2195(17) 


0.2960(35) 


0.7413(88) 


0.05169(24) 


0.916(28) 


6.5 


0.14925 


0.0782(39) 


0.189(13) 


0.414(31) 


0.00860(42) 


0.95(14) 




0.14901 


0.1119(32) 


0.2079(79) 


0.538(21) 


0.01759(42) 


0.951(66) 




0.14877 


0.1394(29) 


0.2232(60) 


0.625(18) 


0.02658(41) 


0.938(47) 




0.14853 


0.1632(25) 


0.2368(45) 


0.689(14) 


0.03565(39) 


0.929(37) 
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TABLE III. B parameters obtained with the WI, Wlys, and PT methods for (3 = 5.9- 6.5 in 
quenched QCD. Operators are renormahzed at 1/a in the NDR scheme. 



B;^(NDR,l/a) 5^(NDR,l/a) 



f3 


K 


WI 


WIvs 


PT 


WI 


PT 


5.9 


0.15862 


0.270(75) 


0.108(26) 


-0.631(13) 


0.0076(21) 


-0.01767(30) 




0.15785 


0.507(39) 


0.259(13) 


-0.0728(56) 


0.0236(18) 


-0.00340(25) 




0.15708 


0.620(26) 


0.3617(82) 


0.1937(38) 


0.0408(17) 


0.01274(26) 




0.15632 


0.687(20) 


0.4348(60) 


0.3501(32) 


0.0585(17) 


0.02977(30) 


6.1 


0.15428 


0.62(14) 


0.226(33) 


-0.416(24) 


0.0156(34) 


-0.01039(49) 




0.15381 


0.686(75) 


0.331(20) 


-0.044(13) 


0.0293(32) 


0.00186(57) 




0.15333 


0.720(51) 


0.406(14) 


0.257(10) 


0.0439(31) 


0.01565(66) 




0.15287 


0.747(39) 


0.461(11) 


0.3813(86) 


0.0585(31) 


0.02985(74) 


6.3 


0.15131 


0.66(16) 


0.286(41) 


-0.157(22) 


0.0175(42) 


-0.00418(53) 




0.15098 


0.713(92) 


0.387(26) 


0.173(12) 


0.0315(41) 


0.00765(56) 




0.15066 


0.745(68) 


0.455(19) 


0.3362(96) 


0.0460(43) 


0.02072(63) 




0.15034 


0.775(54) 


0.511(16) 


0.4449(87) 


0.0625(45) 


0.03586(75) 


6.5 


0.14925 


0.69(39) 


0.212(87) 


-0.338(61) 


0.0115(65) 


-0.00564(87) 




0.14901 


0.65(17) 


0.327(46) 


0.148(25) 


0.0223(59) 


0.00512(92) 




0.14877 


0.67(11) 


0.406(33) 


0.329(20) 


0.0348(58) 


0.0171(12) 




0.14853 


0.699(82) 


0.467(24) 


0.440(15) 


0.0483(57) 


0.0304(12) 



39 



TABLE IV. Results for ^^(NDR, 2GeV) obtained with the WI, WIvs and PT methods as a 
function of (5. Values of S^(NDR,2GeV) in the chiral limit for the WI and PT methods are also 
given. 



13 



5.9 



6.1 



6.3 



6.5 



a = 



5k(NDR, 2GeV) 



5^(NDR, 2GeV) 



mo=0 



WI +0.360(60) +0.66(11) +0.71(12) +0.69(19) 

WIvs +0.162(20) +0.278(27) +0.346(35) +0.347(55) +0.562(66) 

PT -0.391(13) -0.167(20) +0.037(19) +0.180(36) +0.622(69) 

WI -0.0166(32) -0.0055(45) -0.0007(66) +0.0038(96) 

PT -0.03761(67) -0.03055(90) -0.0222(11) -0.0180(15) -0.0023(27) 
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TABLE V. Results for and Ap^(p). In the entries labeled by "color" the color factors 

are given with the use of Cp = 4/3 and Nc = 3. In the third column the label F refers to the 
results in the Feynman gauge A = 1. The label L refers to the coefficients of the gauge parameter 
— (1 — A) for the results in a general gauge. Divergent part is expressed as 1/e = 1/e — 7e + ln47r. 

represents p^7^(l — 75) (8)^1^71^(1 — 75). Definition of evanescent operators E are given in 



(a) NDR 




color 








(a) 


F 


-1/e - Inl^Vp^l - 1 




(a) 


L 


-1/e - Inl^Vp^l - 1 


AS75(p) 


color 




CpSij 




(b) 


F 


7/.75(l/e + In + 1) - '^Piiilb/p^ 




(b) 


L 


7/.75(l/e + In + 1) - '^vd^hlv^ 




color 




CpSij 




(b) 


F 


75(4/e + 4In|^VP^| + 6) 




(b) 


L 


75(l/e+ ln|/iVp'| +2) 


l{p) 


color 




Cp6ij6ki/2 + {5ij5ki - Sii6kj/Nc)/4: 




(c), (d) 


F 


ijl ® ijli^l^ + In |/uVp^| + 1) - 2^^^ ® i^lv^ 




(c), (d) 


L 


l]l ® l\l{^r^ + In \li^lv^\ + 1) - ® i^lv^ 




color 




{6ii6kj - 5ij5ki/Nc)/4: + {5ij6ki - 6ii6kj/Nc)/4 




(e), (f) 


F 


ij: ® 7^'(-4/e - 4 In _ 9 + 3) + E^^P^^Je 




(e), (f) 


L 


7^®7;^'(-l/e- ln|^Vp2| _2 + 4in2) 




color 




- 5ij5ki/Nc)/'i + Cp5ii6kj/2 




(g), (h) 


F 


lli ® 7m + In ImVp'I) + 2^"- ^''/p^ + 




(g), (h) 


L 


ln|/xVP^|) + 2^^«>#^/p^ 
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(b) DRED 



E(l)(p)/(i;f^) 


color 








(a) 


F 


-1/e- ln|/xVp2| _2 




(a) 


L 


-1/e - ln|//2/p^| - 1 




color 








(b) 


F 


7;.75(l/e + In \^lVp'\ + 5/2) - 2p^^^75/p' - ^^/e 




(b) 


L 


7/,75(l/e + In \iJ?/p'^ \ + 1) - 2p^^-i^/p'^ 




color 








(b) 


F 


75(4/e + 41n|/xVp2| +8) 




(b) 


L 


75(l/e+ ln|^Vp2| 


l(X ^ IfJ, 


color 




CFSij6ki/2 + - 6ii6kj/Nc)/4: 




(c), (d) 


F 


7^ ® 7^(1/6 + In I/xVp'I + 5/2) - 2^^^ ® |^^/p' - 

" " 1 fj, ^ 1 fj, 




(c), (d) 


L 


7^ ® 7;^'(l/e + In |/uVp^| + 1) - 2^^ j^^/p^ 




color 




(^^(^ikj - 5ij5ki/Nc)/4. + (5ij(5jk/ - SiiSkj/Nc)/4: 




(e), (f) 


F 


In ® Ini-^/^ - 4 In 1 itVp^l - 8 + 16 In 2) 




(e), (f) 


L 


7^'® 7m - In ^VP^I - 2 + 41n2) 




color 




- ^ij^ki/Nc)/'^ + CFSiidkj/2 




(g), (h) 


F 


if: ® 7^'(l/e + In |/xVp2| + 3/2) + 2^^^- ^ ^^/p2 + 




(g), (h) 


L 


7^'»7;^'(Ve+ ln|/xVP^|) + 2#^(»^^/p^ 



TABLE VI. Scheme dependence of renormalization constants Ayv+AA, ^A, and Ap. 

NDR-DRED NDR-RI 
Avv+AA -3 -14/3 + 8 In 2 

Aa -2/3 

Ai- -4/3 16/3 
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1 1 — I — I I I I I I r 



0.0 - 



-0.2 - 



-0.4 



J I I I L 



— I — I I I I I I 1 1 1 — I I I I I 

mixing coefficients - 

(WI) 



J I I I I 



0.01 



0.10 



1.00 



10.00 



2 2 

p ^ 

Fig.l 
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1 — I I I I I I 1 1 — I — I I I I I 



0.0 - 



mixing coefficients 

(NPR) 

••••5-1 



-0.2 - 
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I 



1 + 
i 



II* 



-0.4 



J I I I L 



J I I I I 



0.01 



0.10 



1.00 



10.00 



2 2 

p ^ 

Fig.2 



0.2 



0.0 - 



-0.2 - 



-0.4 



0.01 



11 — r 



I I I I I r 



t 



J I I I L 



1 — I I I I I I 1 1 1 — I I I I I 

mixing coefficients 
(without first term) 
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1.00 



10.00 
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Fig. 3 (a) 
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